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ABSTRACT 


A linear  triangular  finite  element  formulation  was  used 
to  solve  the  steady  state  two  dimensional  conduction  heat 
transfer  equation.  A FORTRAN  IV  computer  program  of  the 
above  formulation,  employing  double  precision  arithmetic 
and  compact  storage  techniques,  was  applied  to  study  heat 
transfer  in  an  internally  finned  rotating  heat  pipe.  Results 
were  obtained  for  water  in  a copper  and  stainless  steel 
condenser  with  varying  outside  heat  transfer  coefficient, 
rotational  speed,  and  fin  angle  (number  of  fins) . 

Numerical  results  of  the  heat  transfer  rate  in  the  copper 
. condenser  were  shown  to  have  less  than  2%  variance  to  those 

obtained  earlier.  A rotating  heat  pipe  designed  to  operate 
with  a small  value  of  outside  heat  transfer  coefficient 
experiences  no  significant  increase  of  the  heat  transfer  rate 
with  an  increase  in  rotational  speed,  since  the  value  of  the 
outside  thermal  resistance  dominates. 

Heat  transfer  rate  continuously  increases  as  the  fin 
angle  decreases.  However,  the  increase  is  only  slight  when 
the  fin  angle  is  less  than  11  degrees. 
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I.  INTRODUCTION 

A.  THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a closed  container  designed 
to  transfer  a large  amount  of  heat  in  rotating  machinery. 

Its  three  main  component  parts  are:  a cylindrical  evaporator, 
a truncated  cone  condenser,  and  a working  fluid  as  shown 
in  Figure  1. 

At  rotation  above  the  critical  speed  of  a rotating  heat 
pipe,  the  working  fluid  forms  an  annulus  in  the  evaporator, 
and  will  be  vaporized  by  heat  addition  in  it.  The  vapor 
flows  toward  the  condenser,  as  a result  of  a pressure  differ- 
ence, transporting  the  latent  heat  of  evaporation  with  it. 
External  cooling  of  the  condenser  causes  the  vapor  on  the 
inner  wall  to  condense  and  release  its  latent  heat  of  evapor- 
ation. The  centrifugal  force  due  to  the  rotation  has  a 
component  acting  along  the  condenser  wall  that  will  act  to 
drive  the  condensate  back  to  the  evaporator  where  the  cycle 
is  repeated. 

In  a conventional  heat  pipe,  the  force  driving  the  con- 
densate back  to  the  evaporator  is  due  to  capillary  action, 
which  poses  a limit  to  its  operation.  However,  the  driving 
force  in  the  rotating  heat  pipe  can  be  designed  to  operate 
in  any  orientation,  including  in  a zero  gravity  field. 


eat  Out 


Figure  1 Schematic  Drawing  of  a Rotating  Heat  Pipe 


B.  UTILIZATION  OF  AN  INTERNALLY  FINNED  CONDENSER 


The  first  theoretical  investigation  of  the  rotating 
heat  pipe  conducted  at  the  Naval  Postgraduate  School,  was 
performed  by  Ballback  [1]  in  1969.  He  studied  the  limi- 
tations of  performance  imposed  on  the  rotating  heat  pipe 
due  to  various  fluid  dynamic  mechanisms.  Using  existing 
theory  and  experimental  correlations,  he  was  able  to  esti- 
mate the  sonic  limit,  boiling  limit,  entrainment  limit, 
and  the  condensing  limit  of  performance. 

1.  The  Sonic  Limit 

When  increasing  the  heat  flux  in  a rotating  heat 
pipe,  it  is  entirely  conceivable  to  reach  a limiting  flow 
rate  of  the  vapor  brought  on  by  a choked  flow  condition  in 
the  pipe.  This  condition  imposes  a limiting  value  on  the 
amount  of  energy  the  vapor  can  transport , thus  reducing  the 
effectiveness  of  the  heat  pipe.  Since  the  heat  transfer 
method  in  the  rotating  heat  pipe  depends  strongly  on  the 
latent  energy  of  the  working  fluid,  the  rate  of  heat  trans- 
fer in  the  heat  pipe  can  be  expressed  as: 


Qt 


(1-1) 


where  mv  = vapor  mass  flow  rate  in  lbm/hr.  From  the  con- 
tinuity equation, 


m_ 


= P„  U„  A 
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(1-2) 


where 


I 


Uv  = velocity  of  the  vapor  in  ft/sec,  and 

2 

A = cross  sectional  area  for  the  vapor  flow  in  ft  . 

The  heat  transfer  rate  becomes 

°t  ’ hfg  (I-3> 

and  the  vapor  velocity  is  considered  to  be  sonic, 

Uv  = c = /g0kRT  (1-4) 

where 

c = sonic  velocity  in  ft/sec 

2 

gQ  = gravitational  constant  = 32.1739  ft-lbm/lbf-sec 
k = ratio  of  specific  heats 
R = gas  constant  in  ft-lbf/lbm  R,  and 
T = absolute  temperature  in  °R. 

2.  Boiling  Limit 

It  has  been  postulated  by  Kutateladze  [2]  that  the 
transition  from  nucleate  to  film  boiling  is  totally  a 
hydrodynamic  process.  With  this  postulation,  he  determined 
the  theoretical  formula  for  predicting  the  burnout  heat 
flux: 
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***....  ^ r - . UMir  ^ 


(1-5) 


r 


1/4 

Qt  = K/p^  Ab  hfg{°  9(pf  " pv>  } 

= constant  value 

= density  of  the  vapor  in  lbm/ftJ 

2 

= heat  transfer  area  in  the  boiler  in  ft 
= latent  heat  of  vaporization  in  Btu/lbm 
= surface  tension  in  lbf/ft 

2 

= acceleration  of  gravity  in  ft/hr 

3 

= density  of  fluid  in  lbm/ft 

The  experimental  data  obtained  by  Kutateladze  suggested  a 
value  for  K in  the  range  of  0.13  to  0.19. 

3 . Entrainment  Limit 

The  flooding  constraint  in  a wickless  heat  pipe  was 
examined  by  Sakhuja  [3] , who  assumed  that: 

a)  The  vapor  and  liquid  are  exposed  to  each  other 
in  a counterflow  movement. 

b)  The  counterflowing  vapor  tends  to  retard  the 
falling  film  of  liquid  through  shear  stress. 

c)  The  liquid  film  remains  stable  and  smooth,  and 
the  shear  stress  is  usually  small. 

d)  The  buoyancy  forces  resulting  from  the  density 
difference  between  the  vapor  and  liquid  are  the 
means  for  maintaining  the  counterflow. 


where 

K 
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e)  The  buoyancy  forces  are  balanced  by  dissipative 
effects  which  are  proportional  to  the  momentum 


fluxes  of  vapor  and  liquid  streams. 


His  resulting  correlation  for  flooding  is 


A C h£  /gD  (p  ,-p  ) p 
x fg  ’ Hf  Hv  Hv 


[1  + (Pv/Pf ) 7 1 


(1-6) 


where 


= heat  transfer  rate  in  Btu/hr 


= flow  area  in  ft' 


dimensionless  constant,  0.725  for  tube  with 
sharp  edged  flange 


= latent  heat  of  vaporization  in  Btu/lbm 


= acceleration  due  to  gravity  in  ft/hr' 


= inside  diameter  of  heat  pipe  in  ft 


= density  of  the  fluid  in  lbm/ft' 


= density  of  the  vapor  in  lbm/ft' 


Condensing  Limit 


Ballback  [1]  determined  the  condensation  solution 


for  a rotating  heat  pipe  by  modeling  the  condenser  section 


of  a rotating  heat  pipe  as  a rotating  truncated  cone,  from 


which  the  following  condensation  limit  was  developed: 


.2  kt  Of  “2  hfa<WV 


sin  <f> 


} { (R  + Lsin  <J>]  - R 


(1-7) 


where 


Qt  = total  heat  transfer  rate  in  Btu/hr 

kf  = thermal  conductivity  of  the  condensate  film 

r in  Btu/hr-ft-F 

pf  = density  of  fluid  in  lbm/ft3 

to  = angular  velocity  in  1/hr 

hfg  = latent  heat  of  vaporization  in  Btu/lbm 

Ts  = saturation  temperature  in  F 

Tw  = inside  wall  temperature  in  F 

= viscosity  of  fluid  in  lbm/ft-hr 

<f>  = half  cone  angle  in  degrees 

R = minimum  wall  radius  in  ft 

o 

L = length  along  the  wall  of  the  condenser  in  ft 

The  condensing  limit  equation  is  a function  of  the  geometry 
and  speed  of  the  rotating  heat  pipe,  and  the  physical 
properties  of  the  working  fluid. 

Tantrakul  [4]  calculated  these  limitations  on  a 
heat  pipe  with  specific  physical  characteristics  as  shown 
in  Table  1,  with  the  results  shown  in  Figure  2. 

TABLE  1 


Specification  of  a Typical 

Rotating 

Heat  Pipe 

Length 

14.000 

inches 

Minimum  diameter 

2.000 

inches 

Wall  thickness 

. 0.125 

inches 

Internal  half  angle 

1.000 

degree 

Rotating  speed 

2700 

RPM 
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Obviously  from  the  results  in  Figure  2,  the  condensing  limit 
line  is  the  predominant  limitation  for  the  amount  of  heat 
that  can  be  transferred  from  the  heat  pipe.  However,  the 
other  limitations  may  become  important  as  the  heat  pipe 
geometry  and  operating  conditions  are  varied. 

In  order  to  augment  the  heat  transfer  capacity  of 
the  heat  pipe,  the  current  efforts  are  aimed  at  raising  the 
condensing  limit  line  which  may  be  accomplished  by: 

a)  a high  value  of  cone  angle,  to  increase  the 
driving  force. 

b)  some  type  of  promoter  of  dropwise  condensation 
to  increase  the  value  of  the  convective  heat 
transfer  coefficient,  h. 

c)  use  of  an  internally  finned  condenser  to  increase 
the  inner  wall  surface  area  and  the  value  of 

h,  since  the  presence  of  the  fin  will  decrease 
the  condensate  film  thickness. 

By  machining  fins  in  the  inner  wall  of  the  condenser, 
a significant  improvement  of  heat  transfer  will  be  realized. 

C.  ANALYSIS  OF  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 

A disadvantage  of  improving  the  performance  of  a rotating 
heat  pipe  with  fins  machined  in  the  inner  wall  is  that  an 
increase  in  the  wall  thermal  resistance  results  due  to  the 
increase  in  wall  thickness  brought  on  by  the  fins.  By  using 
a material  with  a high  value  of  thermal  conductivity,  this 
effect  can  be  minimized. 
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ternally  Finned  Condenser 
ometry  Showing  Fins , Troughs 
and  Lines  of  Symmetry 


Schafer  [5]  developed  an  analytical  model  for  the  case 
of  fins  with  a triangular  profile  as  shown  in  Figures  3 and 
4.  In  the  analysis  he  assumed  one  dimensional  heat  conduc- 


tion through  the  wall.  Corley  [6]  for  this  same  case 
developed  a two  dimensional  heat  conduction  model  using  a 
Finite  Element  Method,  and  also  assumed  a parabolic  tem- 
perature distribution  along  the  fin  surface.  His  results 
indicated  a significant  improvement  in  heat  transfer  per- 
formance of  about  75%,  above  that  predicted  by  the  one 
dimensional  model  of  Schafer  [5].  However,  Corley  [6] 
cautiously  noted  a probable  error  of  50%  existed  in  the 
calculation  of  the  heat  transfer  at  the  fin  apex,  and  conse- 
quently mentioned  that  there  may  be  a total  heat  transfer 
error  of  as  high  as  15%.  Also,  he  noted  that  his  analysis 
was  not  valid  fro  a stainless  steel  condenser  because  of 
numerical  difficulties.  Tantrakul  [4]  modified  Corley's 
computer  program  by  increasing  the  number  of  finite  elements 
in  order  to  minimize  the  apex  heat  transfer  error.  His 
results  with  this  modification  converged  with  the  results 
of  Corley. 

D.  THESIS  OBJECTIVES 

The  objectives  of  this  thesis  are: 

1.  To  generate  a general  computer  program  for  two 
dimensional  steady  state  Conduction  Heat  Transfer. 

2.  To  use  this  program  to  study  heat  transfer  in  a 
finned  rotating  heat  pipe. 
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II.  FINITE  ELEMENT  SOLUTION 


A.  REVIEW  OF  THE  PREVIOUS  ANALYSIS 

Schafer  [5]  studied  the  one  dimensional  model  heat 
transfer  solution  and  Corley  [6]  studied  the  two  dimensional 
model  for  an  internally  finned  rotating  heat  pipe.  Both 
used  the  same  assumptions  and  boundary  conditions  based 
upon  the  analysis  of  Ballback  [1] , which  are  similar  to 
those  used  in  the  Nusselt  analysis  of  film  condensation  on 
a flat  wall.  The  more  important  of  those  assumptions  are: 

1.  steady  state  operation 

2.  film  condensation,  as  opposed  to  dropwise  condensation 

3.  laminar  condensate  flow  along  both  the  fin  and  the 
trough 

4.  static  balance  of  forces  within  the  condensate 

5.  one  dimensional  conduction  heat  transfer  through 
the  film  thickness  (no  convective  heat  transfer  in 
the  condensate  film) 

6.  no  liquid  - vapor  interfacial  shear  forces 

7.  no  condensate  subcooling 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of 
the  wall  section  (symmetry  conditions) , as  shown 
in  Figures  5 and  6 

9.  saturation  temperature  at  the  fin  apex 

10.  zero  film  thickness  at  the  fin  apex 

11.  negligible  curvature  of  the  wall 
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t * 0.0312S  in. 
b *0.025  in . 

R • - 0.762  in. 
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Figure  6 Internally  Finned  Condenser 
Finite  Element  Geometry  Con- 
sidered in  Analysis  of  Corley  [6] 
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The  complete  development  of  Schafer's  finned  condenser 
heat  transfer  is  presented  in  his  thesis.  His  results  are 
used  as  a comparison  in  this  thesis. 

However,  since  the  same  fluid  flow  and  heat  transfer 
theory  applies  to  the  two  dimensional  conduction  model, 

Corley  [6]  utilized  this  model  and  further  developed  it 
using  the  Finite  Element  Method.  He  considered  an  internally 
finned  geometry  as  shown  in  Figure  6,  and  used  the  fin 
condensate  film  equation  derived  by  Schafer  [5] 


k (T_-T  (z) ) u-  dz 
(z)d<5  (z)  = r b r 


2 2 

oj  (Rq  + x sin  $)  h^  cos  <p  cos  a 


(II-l) 


where 

6 = the  fin  condensate  film  thickness  in  ft- 

kf  = thermal  conductivity  of  working  fluid  in 

Btu/hr-f t-F 

= density  of  working  fluid  in  lbm/ft^ 

= viscosity  of  working  fluid  in  lbm/ft-hr 
a)  = angular  velocity  in  1/hr 

Rq  = minimum  radius  in  ft 

x = distance  along  apex  of  fin  in  ft 

Z = distance  along  the  fin,  from  the  apex  in  ft 

<p  = half  cone  angle  in  degrees 

hfg  = latent  heat  of  vaporization  in  Btu/lbm 


He  assumed  a parabolic  temperature  distribution  along  the 
fin  surface, 


Vz)  = riz  + siz  + fci 


(II-2) 


to  calculate  the  convective  heat  transfer  coefficient  along 
the  fin  surface.  He  developed  a solution  of  the  steady 
state  two  dimensional  conduction  heat  transfer  equation: 


Vk  • 7T  = 0 


(II-3) 


with  boundary  conditions  as  shown  in  Figure  5 
a)  saturation  temperature  at  the  apex. 


T = T 


b)  convective  boundary  condition  in  region  [1] 


-k  il 

3n 


h1(T  - T ),  and  in  region  [3] 


h2 (T  - J 


c)  symmetry  conditions,  = 0,  in  region  [2]  and 
region  [4] . 

He  used  the  8 node  quadratic  isoparametric  element  in  the 
mesh  shown  in  Figure  6,  and  a modification  of  the  three 
dimensional  isoparametric  Finite  Element  Method  computer 
program  assembled  by  Lew  [7] . 

Applying  the  boundary  condition  a) , T = Tsat  at  z = 0 
gives  t^  = Tsat  an<*  e<3uat-j-on  (H“2)  becomes: 


Tw(z)  = rxz  + sLz  + Tsat, 


( II-4 ) 


Tsat  - Vz)  = *rlz  - slz' 


(II-5) 


Upon  substituting  equation  (II-5)  into  equation  (II-2)  , 
and  integrating  the  result  from  z = 0 to  z , gives  an 


equation  for  6(z): 


5(z)  = [- 


3 2 

4 kf  (_ri  V " si  V)uf 


p^  00  (Rq  + x sin  <j>)  h^  cos  cos  a 


( II-6 ) 


Assuming  one  dimensional  conduction  heat  transfer  through 
the  thin  condensate  film,  the  local  convective  heat  transfer 


coefficient  is: 


h (z) 


3 2 2 

k p uj  (Rn  + xsin  cf> ) h ^ cos  $cos  a 1/4 
r r r u rg  i 


, , z z x 

4 -J-) 


(II-7) 


where  kf  = fluid  thermal  conductivity  in  Btu/hr-ft-F  which 
is  applicable  from  z = 0 to  z = Z*  = ZQ  - 6*/cos  a ; where 
6*  is  the  trough  condensate  film  thickness,  as  shown  in 
Figure  5.  Across  the  width  of  the  trough,  the  convective 
heat  transfer  coefficient  is: 


V 


h (x) 


trough 


kf 

VHJY 


(II-8) 


The  heat  transfer  rate  for  the  entire  condenser  is  then 
determined  within  the  computer  program  algorithm  by  dividing 
its  axial  length  into  100  increments  of  length  Ax  (see 
Figure  4) , and  solving  for  the  incremental  heat  transfer 
rate  using  an  iterative  procedure.  The  incremental  rates 
are  then  summed  to  yield  the  total  heat  transfer  rate. 

Figure  7 shows  a flow  chart  description  of  the  computer 
program  used  in  Corley's  analysis  [6],  Constant  values  of 
r^  and  s^  are  determined  at  each  increment  using  Lagrangian 
interpolation  fit  of  the  temperature  at  the  nodes  1,  3, 
and  6 (see  Figure  6)  from  the  previous  iteration.  At  the 
first  increment,  values  of  r^  and  s^  are  calculated  from  an 
initial  guess  of  the  temperature  at  these  nodes.  The  average 
nodal  values  of  h^aVg/  ^3aVg  an<*  the  average  of  h of  the 
left  hand  boundary  region  associated  with  node  6 (see  Figure 
6),  are  calculated  using  the  position  dependent  value 
(equation  (II-7))  in  a Simpson's  Rule  approximation.  In 
the  trough,  h(x)  trough  ^"s  calculate<^  from  equation  (II-8), 
using  the  6*  value  calculated  from  the  previous  increment, 
and  this  value  is  used  as  the  average  convective  heat  trans- 
fer coefficient  for  nodes  9 and  12,  and  the  right  hand 
boundary  of  node  6.  The  average  heat  transfer  coefficients 
calculated  for  the  left-  and  right-hand  boundary  region  of 
node  6 are  weighted  by  their  respective  region  length.  These 
values  are  then  added  to  yield  the  average  convective 
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Figure  7 Computer  Program  Flowchart  of  Two 
Dimensional  Conduction  Analysis  of 
Corley  [6] 
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The  remaining  boundary 
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coefficient  for  node  6,  hgaVg* 
conditions  are  now  applied,  and  the  incremental  heat  trans- 
fer rate  per  unit  of  condenser  length  Qi,  and  the  nodal 
temperature  values  are  generated  from  the  Finite  Element 
Method  subroutine.  If  the  resultant  heat  transfer  rate 
converges  to  within  a 0.05%  change  from  the  value  of  the 
previous  iteration,  i.e. 


< 0.0005 


(II-9) 


then  the  increment  is  deemed  solved  and  Qi  is  set  equal  to 
the  final  iteration,  Q ^ . However,  if  the  convergence 
criterion  is  not  met,  newly  calculated  nodal  temperature 
values  are  used  to  establish  a new  set  of  inner  wall  convec- 
tive heat  transfer  coefficients,  and  the  Finite  Element 
solution  is  repeated  until  the  heat  transfer  rate  converges 
to  within  0.0005  tolerance.  The  incremental  mass  flow  rate 
in  the  trough  is  now  calculated  using  the  relation 


M 


tot  i 


( 11-10 ) 


where  this  value  is  added  to  the  mass  flow  from  the  previous 
increments.  A modification  of  the  trough  mass  flow  rate 
equation  from  Schafer's  thesis  [5]  is  used  to  calculate  the 
subsequent  interval's  trough  condensate  thickness  6*(x), 
with  the  polynomial  root  finder  subroutine: 
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1 1 fan  jfr 


2 2 2 
p w (R-  + xsin  0)6*(.x)  sin  <j>  ? 

(x)  = * [6*  (x)  e+6*  (x)  tan  a] 


tot 


(11-11) 


This  incrementation  is  continued  until  the  entire  length  of 
the  condenser  is  transversed.  Incremental  heat  transfer 
rates  are  summed,  hence  the  total  heat  transfer  rate: 


tot 


100 

2 Nf in  l (QiAx) 
i =1 


(11-12) 


where  is  the  number  of  fins  along  the  inner  wall 

circumference . 

To  start  the  incrementation,  two  initial  guesses  are 
made.  These  are  the  initial  value  of  5*  which  is  taken 
from  the  analysis  by  Sparrow  and  Gregg  for  condensation  on 
a rotating  disk  [8]  and  the  initial  values  of  the  tempera- 
tures at  the  nodes  1 , 3 and  6 . 

Corley  [6]  noted  that  an  error  of  50%  due  to  calculation 
of  heat  transfer  rate  at  the  apex  may  cause  a total  error 
of  15%.  To  minimize  this  error  he  suggested  an  increase  in 
the  number  of  finite  elements. 

Tantrakul  [4]  continued  Corley's  study.  He  increased 
the  number  of  finite  elements  between  the  apex  of  the  fin 
and  the  trough  as  shown  in  Figures  8 and  9 . To  save  com- 
puter execution  time,  he  reduced  the  number  of  increments 
to  25  instead  of  100  increments  as  in  Corley's  calculation. 
However,  even  with  the  reduced  number  of  increments,  the 
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Figure  9 Internally  Finned  Condenser  Finite 

Ele  ent  Geometry  (4  Finite  Elements) 
Considered  in  Analysis  of  Tantrakul  [4] 
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results  of  Tantrakul  [4]  converged  favorably  with  Corley's 
results  [6]  . 


B . FORMULATION 

In  this  thesis  the  author  used  the  same  physical  rela- 
tion and  technique  to  calculate  the  heat  transfer  rate  from 
the  condenser.  These  are  listed  below: 

a)  assume  parabolic  temperature  distribution  along  the 
fin  boundary 

b)  utilize  equation  (II-6)  to  determine  the  condensate 
thickness,  6(z) 

c)  use  equation  (II-9),  (11-10) , and  (11-11)  to 
calculate  the  mass  flow  rate  in  the  trough,  and 
determine  the  subsequent  condensate  thickness. 

d)  iterate  to  arrive  at  the  desired  solution  with  the 
convergence  criterion  0.01,  since  the  results  are 
converged  at  this  value  (see  Table  2)  . 

For  the  Finite  Element  Method  solution,  the  author  used  a 
linear  triangular  finite  element  model  as  shown  in  Figures 
10  through  12. 

Corley  [6]  solved  the  field  equation  for  the  fin  with 
the  following  boundary  conditions: 

a)  temperature  at  the  apex  is  equal  to  the  saturation 
temperature 

b)  convective  boundary  conditions  in  the  inside  and 
outside  of  the  condenser 

c)  symmetry  condition  on  the  sides. 
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The  first  boundary  condition  indicates  no  heat  transfer 
at  the  apex  and  a value  of  infinity  for  the  convective  heat 
transfer  coefficient  h.  Since  the  value  of  h in  reality 
never  equals  infinity,  the  author  let  the  value  of  tempera- 
ture at  the  apex  float,  which  is  in  agreement  with  the 
theoretical  analysis  for  the  extended  surface  in  a one 
dimensional  case.  Only  for  a fin  with  a small  angle  or  a 
very  long  length  can  the  temperature  at  the  apex  become 
equal  to  the  environment  temperature  (saturation  temperature 
in  this  case) . 

Therefore  the  statement  of  the  problem  for  the  formula- 
tion of  the  Finite  Element  Method  is  shown  in  Figure  13: 


32t  + if?  = o 


3x 


(11-13) 


3y 


with  the  boundary  conditions : 


3T 


a)  in  region  1,  -k  = h..  (T  - T . ) (11-14) 

dn  l sat 


b)  in  region  3,  -k  = h2 (T  - TJ 


c)  in  region  2 and  4, 


3T 

■5K 


= 0 


(11-15) 

(11-16) 


Dividing  the  domain  into  linear  triangular  finite  elements, 
the  differential  equation  within  each  element  is: 


3^T^6^  3^T^e^ 

T“ + V 2”  = 0 


(11-17) 


3x 


3y 
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constant 


L^j  h2 

T„ 

3fT  + i!l  a o 

3x2  3y2 

b.c. 

a)  • k Ik  " hl  (T  • Tsat)  in  re8ion  (!) 

b)  ' k 5K  ■ h2  (T  ' T.J  in  roeion  P) 

c)  ■»  0 in  region  [2]  and  f<lj 

Figure  13  Differential  Equation  and  Boundary 
Conditions  Considered  in  Analysis 
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with  the  boundary  conditions : 


a)  -k 


(e)  3T 


h(e)  (T(e)  - T (e)  ) 


(11-18) 


for  elements  with  sides  along  the  convective  boundary,  and 


= 0 


(11-19) 


for  others.  In  the  above  equations. 


= temperature  within  the  element 

= thermal  conductivity  of  the  element;  assumed 
to  be  constant 

= convective  heat  transfer  coefficient  along 
the  boundary  of  the  element,  assumed  to 
be  constant 

= surrounding  temperature  at  the  boundary  of 
the  element. 


Define  the  approximate  value  of  the  temperature  distribution 
within  each  element  as: 


e (x,y)  = l Ni(x,y)  Tj_ 

i=l 


(11-20) 


wnere 


= Two  dimensional  linear  shape  function,  and 
(e) 

T.  = nodal  temperature. 


Using  equation  (11-20) , equation  (11-17)  can  be  written  in 
approximate  form  as 


2 

3 N. 


R = 


l Ti<e>  -ji+  I T. 

3x  i=1 


i=l 


where 

r = residual. 

Applying  the  Galerkin  criterion 

( R N.  dA 

.(e) 


(e)  = 0 


(11-22) 


and  substituting  equation  (11-21)  into  equation  (11-22)  , 
gives : 


/lTi(e)  dAle) 


= 0 


(e) 


(e) 


3y 


Applying  Gauss'  theorem  to  the  factors  on  the  right  side 
yields  the  symmetric  equation: 


where 


The  solution  of 


I 


3N  9N  3N.3N.  , . 

f (- — i__J — + -_±-_l)  dA  e 

> 3x  3x 


(e) 


3y  3y 


(p..)(e)  (11-26) 


where 


b . b . c . c . 

-Vt  + -±-l 

4A 


W 4A(e) 


(11-27) 


bi 


Yi  - Yk 


(11-28) 


c.  = X,  - X.  (11-29) 

i k d 

j ) = 1,  2,  3 

k ' 


as  shown  in  Figure  14.  Recalling  equation  (11-24),  for 
elements  with  an  edge  along  a convective  boundary , 


(e) 


/N.<t>(T(e)  ) dL  (e) 


T. 
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(e)  h 


(e) 


TeT 


(e) 


/N.N. 
‘ i 1 


dL 


(e) 


/ Nj  dL(e)  (11-30) 

Applying  equation  (11-30)  for  elements  with  the  nodal  point 
1 and  2 located  on  the  convective  boundary,  gives 


- T. 
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(11-31) 


where 


pl'p2 


- one  dimensional  linear  shape  function. 


Therefore,  equation  (11-23)  can  be  written: 
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The  element  matrix  equation  is  then  assembled  into  a system 
matrix  equation  with  the  sequencing  based  upon  element  and 
system  nodal  point  connectivity. 

The  heat  transfer  rate  within  each  element  was  calculated 
by 


(e)  T (e) 

Q(e)  = h(«)  ^(e)  Ax(T1.  . Ts(e)} 


(11-32) 


where 


(e ) 

T^  = temperature  at  node  1 

(e) 

T2  = temperature  at  node  2 
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( Q ) 

Ts  = surrounding  temperature 

( 0 ) 

h = convective  heat  transfer  coefficient 

( 0 ) 

= length  between  nodal  points  1 and  2 
as  shown  in  Figure  15. 

C.  THE  COMPUTER  PROGRAM 

The  program  consists  of  a main  program  and  three 
subroutines ; 

a)  the  routine  "COORD"  used  to  define  positions  of 
system  coordinate  points 

b)  the  routine  "FORMAF"  used  to  formulate  the  Finite 
Element  Method  equations 

c)  the  routine  "BANDEC"  as  an  equation  solver  for  a 
symmetric  matrix  which  has  been  transformed  into 
banded  form. 

The  input  data  is  entered  with  the  units  of  inches  for  the 

2 

length,  degrees  for  the  angle,  Btu/hr-ft  -F  for  convective 
heat  transfer  coefficients,  and  in  F for  the  temperatures. 

Input  data  is  arranged  in  6 classifications,  and  is 
shown  in  Table  3 . 
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TABLE  3 


Classification  Of  Input  Data 


Classification 

I 

II 

III 

IV 

V 

VI 


Type  of  input  data 
system  and  element  connectivity 
condenser  geometry 
statement  of  the  problem 
convergence  criterion 
specification  to  determine  the  system 
coordinate  point 

initial  guess  of  the  temperature  at  the 
nodal  points  along  the  fin  and  the  trough 


A detailed  description  of  each  of  these  classifications 
follows : 

Classification  I 

READ  (5,  100)  NEL,  NSNP,  NBAN 


100  FORMAT  (315) 

Explanation:  NEL 

= number 

of  elements 

NSNP 

= number 

of  system  nodal  points 

NBAN 

= system 

band  width 

READ  (5,  200)  (IEL, (ICOR(IEL,I) ,1=1,3) ,IEL=1,NEL) 

200  FORMAT  (415) 

Explanation:  IEL  = the  element  number 

ICOR ( IEL , I ) = system  nodal  point  (NP)  corresponding 

to  NP  I of  element  IEL 
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Note:  - numbering  the  element,  first  for  all  elements  with 

a convective  boundary,  starting  from  the  top  to 
the  bottom,  then  others 

- numbering  the  nodal  point  of  an  element  is 
counterclockwise 

- the  elements  with . convective  boundary  have  nodal 
point  1 and  2 located  on  the  convective  boundary 

Classification  II 

READ  (5,  300)  CL,  CANGL,  RBASE , R2 , THICK,  BFIN 
300  FORMAT  ( 6G10.5) 

Explanation:  CL  = condenser  length  in  inches 

CANGL  = half  cone  angle  in  degrees 
RBASE  = minimum  radius  of  the  condenser  in 


R2 

THICK 


BFIN  = the  height  of  the  fin  in  inches 
READ  (5,  400)  NDIV,  NEST,  NEFB , NBOTI , NBOTF 
400  FORMAT  ( 515) 

Explanation:  NDIV  = number  of  increments  along  the  length 


of  the  condenser 

NEST  = number  of  the  element  on  the  right 
end  of  the  trough 

NEFB  = the  element  number  with  convective 
boundary  located  at  the  base  of  the 
fin 
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% 


NBOTI  = the  element  number  with  convective 
boundary  located  at  the  right  hand 
of  the  bottom  side 

NBOTF  = the  element  number  with  convective 
boundary  located  at  the  left  hand 
of  the  bottom  side 

Classification  III 

READ  (5,  500)  (MRPM ( I ) , 1=1 , 3 ) 

500  FORMAT  ( 315  ) 

Explanation:  MRPM  = number  of  RPM  (w) 

READ  (5,  600)  (TSS (I) , 1=1 , 3) 

600  FORMAT  ( 3G10.5  ) 

Explanation:  TSS  = saturation  temperature  of  the 

working  fluid  (T  . ) 

READ  (5,  700)  TINF 
700  FORMAT  (G10.5) 

Explanation:  TINF  = outside  temperature  (Ta) 

READ  (5,  800)  (HINT ( I) , 1=1 , 3 ) 

800  FORMAT  ( 3G10.5  ) 

Explanation  HINF  = outside  convective  heat  transfer 

coefficient  (ii  .) 

out 

Classification  IV 
READ  (5,  900)  CRIT 
900  FORMAT  ( G10.5) 

Explanation:  CRIT  = convergence  criterion  (0.01) 
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Classification  V 


READ  (5,  1000)  (FANGL(I) , (=1,3) 

1000  FORMAT  ( 3G10.5  ) 

Explanation:  FANGL  = fin  half  angle  (a) 

READ  (5,  1100)  IFF 
1100  FORMAT  ( 15  ) 

Explanation:  IFF  = n-1,  where  n is  the  number  of  rows 

of  the  upper  triangular  fin  section 
READ  (5,  1200)  (KFIN  (I)  ,KFF  (I)  , (-1 , IFF) 

1200  FORMAT  ( 1615  ) 

Explanation  KFIN  = the  number  of  system  nodal  points 

located  on  the  symmetric  boundary 
of  triangular  fin  section,  but  does 
not  include  the  system  nodal  points 
located  at  the  base  of  the  fin  and 
the  apex 

KFF  = the  number  of  system  nodal  points 
located  along  the  fin  convective 
boundary,  but  does  not  include  the 
system  nodal  points  located  at 
the  base  of  the  fin  and  the  apex 
READ  (5,  1300)  DOBF , DOTH,  JTC,  JLC,  JINT,  KT 
1300  FORMAT  ( 2G10.5,  415  ) 

Explanation:  DOBF  = number  of  column  within  the  fin 

DOTH  = number  of  column  within  the  trough 


JTC  = number  of  the  system  nodal  point 
located  at  the  junction  of  the 
symmetry  boundary  and  the  line  of 
intersection  between  the  fin  and 


the  condenser  wall 


JLC  = the  number  of  the  system  nodal 
point  located  at  the  center  of 
system  coordinate 

JINT  = the  numerical  difference  between  the 
two  adjacent  system  nodal  points 
vertically  at  the  condenser  wall 


section 


KT  = the  number  of  rows  within  the  wall 


section 


Note  : as  an  example  for  25  finite  elements  is  shown  in 
Figure  16 


Classification  VI 


READ  (5,  1700)  T(IE) 
READ  (5,  1700)  T(IG) 


1700  FORMAT  (G10.5) 
Explanation: 


T = initial  values  of  the  temperature 
estimate  at  the  nodal  points  1 and 
2,  for  elements  with  convective 
boundary  along  the  fin  and  the 
trough 
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I 


Figure  16  Specification  for  Input  Data  to 
Determine  the  Coordinate  of  the 
System  Nodal  Points 
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D.  TEST  CASES 

Before  the  Finite  Element  Method  program  (routine 
"FORMAT")  was  applied  to  the  condenser  problem,  a series  of 
test  cases  were  analysed.  These  cases  allowed  a comparison 
to  the  results  obtained  from  previous  analyses  and  reference 
9.  These  cases  are  listed  below: 

1.  A rectangular  block  in  a uniformly  convective 
environment  on  three  sides,  with  a specified 
temperature  on  the  non-convective  boundary.  These 
results  are  compared  with  reference  9 as  shown  in 
Figure  17. 

2.  A rectangular  block  in  a uniformly  convective 
environment  on  the  top  and  bottom  side,  insulated 
on  the  other  sides.  A comparison  of  these  results 
with  the  exact  values  and  Corley's  results  [6]  is 
shown  in  Figure  18. 

3.  A rectangular  block  with  triangular  fin  in  a uniformly 
convective  environment  along  the  fin  and  bottom  side, 
insulated  on  the  wall.  A comparison  of  the  results 
with  Corley's  results  is  shown  in  Figure  19. 

The  results  show  good  agreement  with  the  results  of  reference 
9 (case  1) , and  are  comparable  with  the  exact  values  (case  2) . 
Therefore,  "FORMAF"  can  be  used  for  further  applications. 


h = 120  Btu/hr-ft2-°F 


, 29894.  Btu/hr 

Q ® { 

1 (29894.  Btu/hr) 

Note:  numbers  in  parenthesis  refer  to  the  results  of 

reference  9 

Figure  17  Statement  and  Results  of  Test  Case  1 
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h = 1000.0  Btu/hr-ft^-deg  F 


• 7586.  Btu/hr  actual 
Q = { 7552.  Btu/hr 
’ (7527.  Btu/hr) 


e:  Numbers  in  parentheses  refer  to  results  of 

Corley  [6]  . 


Figure  18.  Statement  and  Results  of  Test  Case  2 


III.  NUMERICAL  RESULTS  AND  DISCUSSION 


The  analysis  concerns  itself  with  the  geometry  of  both 
a copper  and  stainless  steel  condenser  as  tabulated  in 
Table  4,  using  water  as  the  working  fluid. 

TABLE  4 

Specification  of  a Typical  Internally 
Finned  Rotating  Heat  Pipe 


Condenser 

length 

= 

9.000 

inches 

Minimum  radius 

= 

0.775 

inches 

Wall  thickness 

= 

0.03125 

inches 

Height  of 

fin 

= 

0.025 

inches 

Cone  half 

angle 

= 

1.0 

degree 

I 


Results  are  given  in  Tables  5 to  18  and  Figures  20  to  28. 

It  can  be  seen  from  the  data  that  the  solution  of  the  Finite 
Element  Method  converges.  This  can  be  seen  by  examining 
data  of  temperature  distribution  within  the  fins  of  25, 

40,  and  108  linear  triangular  finite  elements  (Tables  7 to 
9),  for  the  same  location  (nodes  10,  15,  and  28  respectively). 

The  heat  transfer  rate  solution  is  generally  affected 
in  the  same  manner  as  the  analysis  of  Schafer  [5]  when  the 
outside  convective  heat  transfer  coefficient  (h  t)  > rotational 
speed  (u>)  , fin  angle  (2a)  , and  thermal  conductivity  (k 
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TABLE  9 

Theoretical  Heat  Transfer  Rate  on  a Stainless  Steel  Finned 

Condenser  at  1000  RPM,  for  h . = 1000.  Btu/hr-f t2-deg.F 

out 


No. 

T 

sat 

Q(Btu/hr) 

o r 

fins 

(°F) 

25  EL 

40  EL 

108  EL 

10 

276 

100 

8948.2 

8944.1 

8893.9 

23764.0 

23727.0 

23607.0 

200 

38557.0 

38492.0 

38286.0 

30 

84 

100 

8515.5 

8494.4 

8461.8 

150 

22511.0 

22445.0 

22354.0 

200  . 

36467.0 

36355.0 

36204.0 

50 

40 

100 

8186.3 

8164.6 

8141.6 

150 

21574.0 

21508.0 

21442.0 

200 

34916.0 

34805.0 

34693.0 

TABLE  10 

Theoretical  Heat  Transfer  Rate  on  a Stainless  Steel  Finned 
Condenser  at  3000  RPM,  for  hQut  = 1000.  Btu/hr-ft^-deg.F 


No. 

of 

fins 

T . 
sat 

(°F) 

Q(Btu/hr) 

25  EL 

40  EL 

108  EL 

10 

276 

100. 

9096.4 

9100.4 

9058.1 

24193.0 

24200.0 

24098.0 

39277.0 

39284.0 

39123.0 

30 

■ 34 

100.- 

8818.4  ' 

8810.3 

8786.0 

150. 

23389.0 

23376.0 

23301.0 

200. 

37934.0 

37900.0 

37781.0 

50. 

40 

100. 

8584.1 

8580.7 

8561.1 

150. 

22720.0 

22690.0 

22642.0 

200. 

36823.0 

36768.0 

36689.0 
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TABLE  11 


Theoretical  Heat  Transfer  Rate  on  a Stainless  Steel  Finned 
Condenser  at  1000  RPM,  for  h . = 5000.  Btu/hr-f t^-deg .F 


No. 

of 

fins 

T 

sat 

(°F) 

Q(Btu/hr) 

25  EL 

40  EL 

108  EL 

10 

276 

100. 

26899.0 

26775.0 

26308.0 

70520.0 

68673.0 

si 

113950.0 

113300.0 

110780.0 

30 

84 

100. 

22742.0 

22540.0 

22286.0 

150. 

58879.0 

58288.0 

57582.0 

200. 

94734.0 

93748.0 

92581.0 

50 

40 

100. 

20426.0 

20253.0 

20086.0 

150. 

52480.0. 

51978.0 

51513.0 

200. 

84275.0 

83443.0 

82676.0 

! 


TABLE  12 


Theoretical  Heat  Transfer  Rate  on  a Stainless  Steel  Finned 
Condenser  at  3000  RPM,  for  hQut  = 5000.  Btu/hr-f t2-deg .F 


Fin 

half 

angle 

No. 

of 

fins 

T _ 
sat 

(°F) 

Q(Btu/hr) 

25  EL 

40  EL 

108  EL 

10 

276 

100 

28705.0 

28701.0 

28263.0 

ISO 

• 75738.0 

75734.0 

74492.0 

200 

122610.0 

122620.0 

120570.0 

30 

34 

100 

25604.0 

25534.0 

25300.0 

ISO 

66822.0 

66734.0 

65951.0 

. 

200 

107830.0 

107800.0 

106260.0 

SO 

40 

100 

23673.0 

23555.0 

23410.0 

150 

61488.0 

61127.0 

60726.0 

200 

99069.0 

98464.0 

97803.0 

TABLE  13 


r 


Temperature  Distribution  Within  the  Fin  at  the 

10t^1  Increment. 

(Copper  Condenser  with  25  Finite  Elements 


TABLE  14 


Temperature  Distribution  Within  the  Fin  at  the 

10fch  Increment 

(Copper  Condenser  with  40  Finite  Elements) 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.01 

15 

98.22 

29 

97.90 

2 

98.84 

16 

98.04 

30 

97.82 

3 

98.83 

17 

97.99 

31 

97.79 

4 

98.69 

18 

98.17 

5 

98.68 

19 

98.16 

6 

98.64 

20 

98.14 

7 

98.55 

21 

98.11 

8 

98.54 

22 

98.06 

9 

98.50 

23 

97.96 

10 

98.45 

24 

97.92 

11 

98.42 

25 

97.99 

12 

98.41 

26 

97.98 

13 

98.38 

27 

97.96 

14 

98.32 

2 8 

97.94 

a = 50  deg 

T 

= 70 

°F 

^out 

* 1000  Btu/hr-f t‘ 

2-deg  F RPM  = 1000 

T 

sat 

* 100  °F 

Q 

= 9687.2  Btu/hr 

1 


TABLE  15 

Temperature  Distribution  Within  the  Fin  at  the 

10th  increment 

(Copper  Condenser  with  108  Finite  Elements 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.04 

24 

98.39 

50 

97.90 

2 

98.91 

26 

98.33 

51 

98.05 

3 

98.90 

28 

98.20 

52 

97.96 

4 

98.79 

29 

98.05 

58 

97.85 

6 

98.76 

30 

97.99 

60 

97.84 

7 

98.68 

31 

97.96 

61 

97.97 

8 

98.68 

32 

98.27 

62 

97.96 

10 

98.63 

34 

98.25 

63 

97.95 

11 

98.58 

36 

98.20 

64 

97.94 

13 

98.56 

38 

98.12 

65 

97.92 

15 

98.51 

39 

98.02 

66 

97.90 

16 

98.49 

40 

97.96 

67 

97.90 

18 

98.47 

41 

97.94 

68 

97.88 

19 

98.45 

42 

98.14 

69 

97.82 

21 

98.37 

44 

98.04 

70 

97.78 

22 

98.41 

48 

97.91 

71 

97.76 

a = 50  deg 

T = 70° 

F 

^out 

= 1000  Btu/hr-ft2 

-deg  F 

RP.M  = 1000 

T . 
sat 

= 100  °F 

Q 

= 9678 

. 9 Btu/hr 

TABLE  16 


Temperature  Distribution  Within  the  Fin  at  the 

10^  Increment 

(Stainless  Steel  Condenser  with  25  Finite  Elements) 


NP 

T ( °F) 

NP 

T(°F) 

1 

99.79 

15 

95.17 

2 

99.08 

16 

94.44 

3 

99.20 

17 

93.94 

4 

99.19 

18 

93.91 

5 

98.24 

19 

93.83 

6 

98.37 

20 

93.69 

7 

97.23 

21 

93.11 

8 

97.21 

9 

97.10 

10 

96.66 

11 

95.48 

12 

95.50 

13 

95.47 

14 

95.36 

a = 50  deg  T = 70°F 

2 

hQut  = 1000  But/hr-f t -deg  F RPM  = 1000 
Tsat  = 100°F • Q = 8186.3  Btu/hr 


TABLE  17 

Temperature  Distribution  Within  the  Fin  at  the 

10t'1  Increment 

(Stainless  Steel  Condenser  with  40  Finite  Elements) 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.82 

15 

96.70 

29 

93.62 

2 

99.31 

16 

95.73 

30 

93.26 

3 

99.39 

17 

95.50 

31 

93.12 

4 

98.67 

18 

95.50 

5 

98.72 

19 

95.48 

6 

98.84 

20 

95.41 

7 

97.97 

21 

95.30 

8 

97.99 

22 

95.12 

9 

98.03 

23 

94.62 

10 

98.09 

24 

94.45 

11 

97.25 

25 

93.93 

12 

97.24 

26 

93.91 

13 

97.20 

27 

93.85 

14 

97.08 

28 

93.75 

a - 50 

deg 

T 

= 7 0 

’F 

^out 

1000  Btu/hr-ft2- 

■deg  F RPM  * 1000 

T = 

sat 

100°F 

Q 

= 8164.6  Btu/hr 
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TABLE  18 


Temperature  Distribution  Within  the  Fin  at  the 

10th  Increment 

(Stainless  Steel  Condenser  with  108  Finite  Elements) 


NP 

T(°F) 

NP 

T(°F) 

NP 

T(°F) 

1 

99.84 

24 

97.21 

50 

91.46 

2 

99.50 

26 

97.13 

51 

94.39 

3 

99.54 

28 

96.64 

52 

94.64 

4 

99.09 

29 

95.83 

58 

94.29 

6 

99.19 

30 

95.52 

60 

93.83 

7 

98.65 

31 

95.43 

61 

93.76 

8 

98.67 

32 

96.32 

62 

93.87 

10 

98.81 

34 

96.28 

63 

93.84 

11 

98.19 

36 

96.15 

64 

93.79 

13 

98.23 

38 

95.82 

65 

93.73 

15 

98.36 

39 

95.31 

66 

93.64 

16 

97.70 

40 

95.03 

67 

93.55 

18 

97.72 

41 

94.95 

68 

93.31 

19 

97.73 

42 

95.45 

69 

93.31 

21 

97.73 

44 

95.41 

70 

93.13 

22 

97.22 

48 

95.04 

71 

93.07 

a = 50 

deg 

T 

= 70  °F 

^out 

1000  Btu/hr-ft^- 

-deg  F RPM  = 1000 

T = 

— sa± 

100°F 

Q 

= 8141 

6 Btu/hr 

Q (Btu/hr) 


1.  h . = 5000  Btu/hr-ft2- deg  F 

out 

2.  h = 1500  Btu/hr-ft2- deg  F 

out 

3.  h « 1000  Btu/hr-ft2- deg  F 

out 

Figure  20  Heat  Transfer  Rate  (Q)  of  Internally 

Finned  Condenser  at  a Particular  Value 
of  h , vs.  Saturation  Temperature  of 
out  Working  Fluid 
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Saturation  Temperature  (°F) 


. a = 30< 


. a = 50° 

Figure  21  Heat  Transfer  Rate  (Q)  of  Internally  Finned 
Condenser  at  a Particular  Value  of  Fin  Half 
Angle  vs.  Saturation  Temperature  of  Working  Fluid 


Q (Btu/hr) 


) 2500  5000  7500  10000  12500  -.15000  17500 

RPM 

hQut  = 1000.  Btu/hr-f 1 2-degF  4.  hout  = 13000.  Btu/hr-f t2-degF 

h . = 5000.  Btu/hr-f t2-degF  5.  h . = 17000.  Btu/hr-f t2-degF 

out  out 

hQut  = 9000.  Btu/hr-f t2-degF 

Figure  23  Heat  Transfer  Rate  (Q)  of  Internally  Finned 

Stainless  Steel  Condenser  vs.  RPM 
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Figure  24  Thermal  Resistance  of 
Copper  Condenser 

Internally 
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Finned 
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Figure  25  Thermal  Resistance  of  Internally  Finned 
Stainless  Steel  Condenser  vs.  h . 
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Q0  : heat  transfer  rate  of 
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smooth  condenser 


Figure  27  Comparison  of  Heat  Transfer  Rate  (Qp/Qg) 
r at  h ^ = 1000.  Btu/hr-ft2-degF 


vs . 


were  varied.  In  Figure  20,  the  influence  of  the  outside 
convective  heat  transfer  coefficient,  h . , is  shown  both 
for  copper  and  stainless  steel.  Figure  21  illustrates  how 
the  rate  of  heat  transfer  was  affected  when  the  fin  half 
angle  was  varied.  In  addition,  the  heat  transfer  rate  was 
found  not  to  vary  significantly  with  rotational  speed  when 
the  heat  pipe  operated  with  low  values  of  outside  convective 
heat  transfer  coefficient,  hQut*  However,  for  high  values 
of  hQut#  the  effect  of  rotational  speed  cannot  be  ignored 
(see  Figures  22  and  23) . 

The  small  variance  of  the  heat  transfer  rate  with  rota- 
tional speed  is  due  to  the  dominancy  of  the  outside  thermal 
resistance,  Rout/  such  that  changes  in  the  condensate  film 
thickness,  6,  will  not  significantly  affect  the  total  thermal 
resistance  as  shown  in  Figures  24  and  25.  The  heat  transfer 
rate  will  have  a large  variance  when  the  heat  pipe  is  operated 

with  a high  value  of  h 

^ out 

From  the  results  of  the  parametric  study,  as  shown  in 
Figures  27  and  28,  the  heat  transfer  rate  continuously  in- 
creases as  the  fin  half  angle  decreases  (b/a  increases) . 

Since  the  total  number  of  the  fins  in  the  condenser  is 
larger  with  a small  half  angle  fin,  this  results  in  an  in- 
crease in  the  heat  transfer  rate.  However,  the  increase 
is  only  slight  when  the  fin  half  angle  less  than  11  degrees 
(b/a  greater  than  5) . 

The  study  shows  the  enhancement  in  heat  transfer  is 
larger  when  operating  with  a high  value  of  hQut  or  when 
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is  small.  Therefore  the  design  operating  with  a 


Rout 

high  value  of  h t#  using  both  internal  and  external  fins 
is  recommended. 
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IV.  CONCLUSIONS 


1.  The  Finite  Element  Method  works  and  converges. 

2.  The  designer  has  a choice  of: 

- material 

- number  of  fins  and  fin  angle 

- RPM 

- h 

out 

3.  The  heat  transfer  rate  increases  continuously 
with  decreasing  fin  angle;  however,  for  half 
angle  less  than  11  degrees,  the  increase  is 
only  slight. 


4 


i 


V.  RECOMMENDAT IONS 

1.  Manufacture  internally  finned  heat  pipe  and  test. 

2.  Manufacture  internally  and  externally  finned  heat 
pipe  and  test. 
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APPENDIX 


COMPUTER  PROGRAM  LISTING 


XXtXXXXXrtXX.XXXXXXXXXXXXXXtXXtXXXXXXtXtXtXXXXXXXXXXXXXXXXXXXXX 


* 

* * 

i TWO  DIMENSIONAL  HEAT  CONDUCTION  FINITE  ELEMENT  ANALYSIS  * 

* OF  ROTATING  HEAT  PIPE  . USING  TRIANGULAR  ELEMENT  MODEL  * 

X COMPILED  BY  MAJOR  IGNATIUS  S PURNOMO  IN  JUNE  ,1973  * 

♦ * 

* * 


xxxxxxii  txxxtxxv.txxxx'irxxtxtxxxxxxxxtxxxxxxxxxxxxxxxxxxxxxxxxxx 

IMPLICIT  REALtSCA-H.O-Z) 

DIMENSION  Z C 200  ) , EPS  ( ECU  ) , HZ  ( 200 ) . XCOF  ( 5 ) . COF 1 5 ) . ROOTr  ( 4 ) , ROOTI  ( * 

> 

* , Tf  200  > • OB  c 200 1 ■ DI1D0T  < 200  ) , UF  ( 200  ) , CF  ( 200  ) , RHOF  ( 200  ) , DEL  ( 200 ) , CU* 

2 

*00 ) ■ AMTOT ( 200  ) , R f 200  ) • Q INC ( 200 ) , TB < 200 ) , TTC 200 ) , TBM ( 200 ) , TEC 200  )~ 

T 

*IBf 200  ) .NCC2001 

DIMENSION  TS3C3 > , FANGLC3  ) , HINFC3 ) , MRPMC  3 ) 

COMMON/’ APOL'yDOBF , DOTH , KF INC  50  ) , K.FF  ( 50  ) , IFF . JTC , JLC , J INT , KT 
COMMON/MaFO/A  C 200  - 50  ) . F ( 200 , 1 ) , H C 200  ) , TS  C 200  ) . TSAT  - CK , NEL , NSNP , N~ 
4 

*N • I COR  C 200 • 3 ) 

COMMON/'PCRD-'X  C 200  ) , V < 200  ) . EZEPO - BU IN . THICK , TALFA , APS 
HDENCAl,Bl.ZZ)-C-1.0DO*CAl*ZZ**3-/3  0D0+Bl*ZZ*X2/2  0D0  n 

ELEMENT  CONNECT IU IT I ES 

READ C 5, 100)  NEL. NSNP. NBAN 
100  FORMAT C 3 IS) 

URITECS.1S0)  NEL. NSNP. NBAN 

ISO  FORMAT C/2X. 'NO  OF. ELEMENTS- '. 15. 10X. 'NO  . OF  SYSTEM  N . P . - ' . 15. 10X . * 

N 

*0  OF  BANDED-'. 15) 


-■  ^ 


ooooooo  o o o o o o o 


200 


READ'  S. 200  ) < IEL, (ICORf IEL, I >, I- 1 , 3 ) . IEL-1 , NEL ) 

FORMAT  < 415 ) 

URITE1'. 6 . 250  ) 

250  FORMAT < /2X  • ELEMENT  , 1CX.  'NP1',14X,  'NP2',15X.  'NP3'  ) 

WRITE*  6. 251  KIEL,  (ICOR(IEL.  I >,I-1.3>.  IEL-i.NEL) 

251  FORMAT ( 15. 12X. IS. 12X . 15 . 12X . 15  ) 

THE  CONDENSER  GEOMETRY 

CL  IS  THE  CONDENSER  LENGTH 

CANGL  IS  THE  CONE  HALF  ANGLE 

RBASE  IS  THE  CONDENSER  RADIUS  AT  THE  BASE 

THICK  IS  THE  UALL  THICKNESS  OF  THE  CONDENSER 

READ < 5 • 300  ) CL . CANGL . RBASE , R2 , THICK . BF IN 
300  FORMAT fbG 10  5) 

CL-CL'12  ODO 
R2*R2,'12  ODO 
RBASE-RBASE/12  ODO 
B*  'IN-BF IN/12  ODO 

NDIU  IS  THE  NUMBER  OF  INCREMENT 

NEST  IS  THE  ELEMENT  NUMBER  AT  THE  END  OF  THE  TROUGH 
NEFB  IS  THE  ELEMENT  NUMBER  AT  THE  BASE  OF  THE  FIN 
NBOTI  IS  THE  FIRST  ELEMENT  AT  THE  BOTTOM  SIDE 
NBOTF  IS  THE  LAST  ELEMENT  AT  THE  BOTTOM  SIDE 

READ( 5 ■ 400  ) NDIU . NEST, NEFB • NBOTI , NBOTF 
400  FORMAT (5 15) 

DIU-DFLOAT(NDIU  ) 

PI-3  1415926335897900 
PHI-2  0D0*CANGL*PI/360  ODO 
SPHI-DSINCPHI  ) 


ooo  ooo  ooo 


DELX-CL^BIU 

CEASE *2 . 0D0*PI CREASE 

REXIT-PBASE+CLTTPHI 

CEXIT-2  ODOtPUPEXIT 

THICK-THICK/12.0D0 

URITE  < 6 , 450 )CL , CANGL . R2 . THICK . RBrtSE - BUIN 
450  FORMAT*.'  ' 'CONDENSER^S  LENGTH- ' . E12  5 . 30X,  'HALF  CONE  ANGLE-' 

m ’ 

* . E12  5, s-SX. 'R2- ' » E12 . 5. 45X. 'THICK- ' . E12  5. / . 5X . 'R-BASE- ' . E12 . 5 . * 
41 

*X. 'B-FIN«',E12  5) 

DATA  FOR  RUNNING 

READ  ( 5 500)  (MRPMC I ) . 1-1 . 3 ) 

500  FORMAT' 315) 

READ (5. 600  ) (TSS(I). 1-1.3) 
eeo  FORMATf  3G10  5) 

PEpQc 5 .700)  TINF 
700  FORMAT fG10  5) 

READ (5. 800  > (HINF(I). 1-1,3) 

800  FORMAT  C 3G10 . 5 ) 

THE  CONUERGENCE  CRITERIAN 

READ (5. 900)  CRIT 
900  FORMAT (G10  9) 


INTERNALLV  FIN  GEOMETRY 


nooo 


READ'S,  1200)  <KFIH<  I ) , ICFF ( I ) » I«1 . IFF  ) 

1200  FORMAT ( 1615 ) 

READ*.  5 - 1300  >DOBF  , DOTH . JTC , JLC , JINT , KT 
1300  FORMAT i £010. 5. 415) 

MHB*NEFB/'2 

NTM  » HBOT I + ( NBOTF-NBOT I )/2 
MtF-HBOTF+1 

WRITE ( 6 . 1350)  ICORCNBOTI .2), ICOR<NEFB. 1 ) . I COR (NTM. 2 ) . ICORCNEST , 1~ 

>, 

•TICOP'NBOTF  1 ) 

1350  FORMAT  ( ///5X , 'TIB* ' . 15, 10X,  'TT*  M5, /.5X,  'TBM« ' > 15, 10X.  'TE«'.I5.~ 

*6X.'TB-' .15) 

DO  1400  KIA-1.3 

ALFA-FaN6L(KIA >*2  OD0*PI/36O  0DO 
SALFA-DSIWALFA ) 

CALFA»DCOS< ALFA ) 

TALFh-DTaNi ALFA) 

EZERO-2  ODO*BUIN*TALFA 
EPSO-EZERO 

NFIN-CBASE.'(EZEPO+EPSO  ) 

EPSEX- rCEXIT-(NFIN*EZERO ) )/NFIN 
BETA- ( EP5EX-EPS0  >/DIU 
ZZERO-BOIN/CALFh 
ZA-0  0D0 

BOUNDARY  CONDITIONS  AND  TEMPERATURE'S  ESTIMATE 
ALONG  THE  FIN  BOUNDARY 

DO  1450  NTINF -NBOTI , NBOTF 
1450  TSCNTINF )«TINF 

DO  1500  NNT-N3F.NEL 
TSCNNT )*0  0D0 
1500  H(NNT)-0  0D0 
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ooo 


DO  1600  IGT« 1 , NEST 
IE-IC0R'IGT.2> 

1600  READ'.  S-  1700  ) T<!IE) 
IG-ICORCNEST, 1 ) 

PEHt'CS.  1700)  TCIG) 

1700  FORMAT* G10  5) 

DO  HOO  KIB-1,3 
MRPM*riPPrifKIB ) 

0MEGA-NRPM*2  ODOtPI*60  ODO 
DO  1400  KIC-1,3 
DO  1200  KL-NBOTI.NBOTF 
1300  HO-  L)*HiriF(KIC  ) 

HIFN-HINFCKIC 1 
DO  1400  KID-1.3 
TSaT-TSSCKID ) 

DO  1900  NSAT-l.NEST 
1900  TS(NSAT  >«TSAT 

TSOLID- (TSAT+TINF )/2  ODO 
DELlD-0  OOOOS7S2D0 
QT-0  ODO 
QTOT-O  ODO 
DMTOT-O  ODO 
NK-NDIU+1 
DO  2000  NI-l.NK 
RCNI >-R2+NI*DELXXSPHI 
EPSCNI  )-EPSO+Nl*B£TA 
APS-EPS(NI  ) 

NODAL  POINTS  COORDINATE 

CALL  COORD 
2(1 )-2A 

DO  2100  IZEL-l.NEFB 
NA-ICORCIZEL. 1 ) 
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ooo  oooo 


NB-IC0P-:iZEL.2) 

XE-tf<HA>-X(MB) 

YE«V«NA)-Y(NB ) 

ELZ-DSGPT*  XE**2+YE**2 ' 

2100  Zv  IZCL+1  )-Z<  IZED+ELZ 

XZB»X< ICORf NHB • 1 ) )-X(ICOR< 1.2)) 

YZE-Vi I COR  < NHB.  1 > )-Y( ICOR( 1.2)) 

ZE«DSQFT<  XZB**2+YZB**2 ) 

ZC-ZZEPO 

1 IPI-1 

PARABOLIC  TEMPERATURE  DISTRIBUTION  ALONG  THE  FIN 
BOUNDARY. USING  LAGRANGE  INTERPOLATION 

2 TPl-TtICORf 1.2 > ) 

TP2«T< I COR (HHB. 1 ) ) 

TP3-Tf ICORCHEFB. 1 ) ) 

AP1-TP1'(ZB*ZC) 

AP2“TP2/'(ZB*(ZB-ZC) ) 

AP3-TP3/(ZC*(ZC-ZB) ) 

BP1— (ZB+ZCUAP1 
BPS— ZC*ftP2 
BP3— ZB*AP3 
A1»APH-AP2-*^.P3 
B1-BPH-BP2+BP3 
TC-0  0D0 

DO  2200  NY -1. NEST 
!00  TC»TC+T(  I COR  (NY,  2)  ) 

AY ■ DFLOAT ( NY+ 1 ) 

TF ■ ( TC+T ( ICCRC NY. 1 ) )+AV*TS(NY ) >/<2  ODOtAY ) 

SOLID-FLUID  PROPERTIES 
HFG-1097  2DO-0  60187SDO*TS<1 > 
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• J . X . . . 


oooo 


RHOF'NI  >»62.774B0-0  0O2S569SD0*TF-O . 0000S3572D0*TF**2 
CF'.III  )*0 . 303400+0  OO0738927PO*TF-0  000001 47321D0*TF**2 
UFCMI )-0 . 001397D0-0  000O14669D0tTF+O  0OO00O0631253D0*TFX*2-0 . 000^ 

*0OO0097G569D0*TFX*3 
UFCNI  )-3600tUFUU  ) 

CUCHI ) =231 . 7772D0-0 . 02222D0XTS0LID 
CK-CU(NI ) 

CONST  * RHOF  ( N I m2*0MEGA**2*HFGXCPHl*CALFAXRCNI  ) 

AUERAGE  ELEMENT  CONUECTIUE  COEFFICIENT  ALONG  THE 
FIN  BOUNDARY 

ZSTAR*ZZERO-DEL(NI VCALFA 
AZZ-DELCNI 5/SALFA 
ZZ=ZSThR 

AZS-DABS( 4*CF ( NI  )*UF (NI  )*HDEN( A1 , B1 . ZZ ) /CONST )**0  25D0 
HAC=0  0DO 
DO  2300  IEL-l.NEFB 
AZ-Z(IEL) 

BZ-ZCIEL+l) 

IF ( ZSTAR  LE  BZ ) GO  TO  2400 
GO  TO  2500 

2400  I F ( H AC  NE  O OD0)  GO  TO  2600 
BZ-ZSTAR 

2500  IFCIEL  NE1)  GO  TO  2700 
AK*  ( BZ-AZ  VS  . 0DO 
ZZ-AK 

GO  TO  2800 

2700  AK -( BZ-AZ J/4.0D0 
ZZ-AZ 


2800  ZEL-4XAK 

DO  2900  NH-1.5 

HZ(NH)-DABS(CF(NI  )**3*C0NST/< 4*UF(NI )*HDEfH A1 . Bl . ZZ ) ) )**0  25D0 


o o o o o o 


2300  ZZ-ZZ+AK 

CONH-AK*  (HZ(  1 )+4»H:(2  )+2*HZ(3  )+4*HZ(4  HHZCS  ) )/<3*ZEL) 
IFiZSTAR  EQ.BZ)  GO  TO  3000 
H(IEL)-CONH 
GO  TO  2300 
3000  AZ-ZSTAR 

HAZ-CONH* i AZ-Z( IZL ) ) 

DELA-AZS 

3100  BZ-Z(IELM) 

DELB- ( BZ-ZSTAR  )*AZZ/ (ZZERO-ZSTAR ) 

DELZ- 1.  DELA+DELB  '/2  ODO 
HAC-(BZ-AZ)*CFCNI >/DELZ 
H(IEL)-(HAZ+HAC)/(BZ-Z(IEL) ) 

GO  TO  2300 
2600  AZ-ZCIEL> 

DELA-DELB 
HAZ-0  ODO 
GO  TO  3100 
2300  CONTINUE 

NETI-NEFB+1 
DO  3200  IEL-NETI . NEST 
3200  HC IEL  )-CF(Nl  i/'DEL(NI ) 

ENTRY  INTO  THE  FINITE  ELEMENT  SOLUTION 

CALL  FORMAF 

CALL  BANDECCNSNP.NBAN. 1 ) 

THE  TEMPERATURE  DISTRIBUTION 

DO  3300  NT-l.NSNP 
3300  T(NT)-F(NT. 1 ) 

TIB(NI  )-T( ICORCNBOTI . 2 ) ) 

TT(NI  )-T( ICOR(NEFB . 1 ) ) 
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TBMCNI  ) «T( ICOR ( NTM .2)) 

TE(NI  >*T ( ICOR' flEST. 1 ) ) 

TE  < HI  • -T ( ICOR  C NBOTF , 1 ) ) 

TTS-0  OD0 
DO  3400  NS-l.NSNP 
3400  TTS-TTS+TCNS) 

PN-DFLOAT'.NS » 

TSOLID-TTS/PN 

C 

C Q AT  THE  BOTTOM  SIDE 

C 

OBI "0  ODO 

DO  3500  IBEL-NBOTI. NBOTF 
NKA-ICORf IEEL. 1 > 

NKB ■ ICOR ( I BEL . 2 > 

XB-XCNKH  i-X(NKB ) 

YB-YCNKA  J-YCNKB ) 

ELB - DSORT  C XBX*2+YB* t2 ) 

3500  QBI-OE I+CT(NKA  )+T<NKB  )-2*TS< IBEL ) )*ELB*H( IBEL )^2 . 0D0 
OB(NI )-QBI*DELX 
C 

c ITERATION  UNTIL  CON'JERGENCE  CRITERIA  IS  MET 


IFCin.EO  1)  GO  TO  3600 

QJ-QBI 

GO  TO  3700 

3600  QI-OBI 
IM-2 
GO  TO  2 

3700  AO-DABS(QJ-QI )/0J 

IF(AQ  LE.CRIT)  GO  TO  3800 

QI-QJ 

GO  TO  2 

3800  DMDOT(NI  )-2*QBI*DELX/HFG 
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ooo 


1 

DMTOT-DHTOT+DnDOTi.NI  ) 

C1*RHQF (HI  )**2*0MEGA**2*R(NI  )/(3*UF(NI ) ) 

XC0F<1  >— DP1T0T 
XCOFt'2  '"0 . ODO 
XCOF(3)-0  0D0 
XC0F<4  >-Cl*EPS».NI ) 

XC0Fi5>-Cl*TALFA 
M-4 

CALL  DPOLPT ( XCOF . COF . M , ROOTR . ROOT I . IER  ) 

IF ( ROOTR (l).LT  0 005D0  AND  RCCTR ( 1 ) GT  0 . 0D0 ) GO  TO  3300 

IF ( ROOTR f 2 ) . LT  O OOSDO  AND . ROOTR ( 2 ) GT  0.0DO)  GO  TO  4000 

IF (ROOTR (3 ) LT  0 C05DO  - .ND  ROOTR ( 3 ) GT  .0  OD0 ) GO  TO  4100 

IF ( ROOTR ' 4 ' LT  0 005D0 . AND  R00TR(4)  GT.0.0DO)  GO  TO  4200 

URITEf  to • 4300  ) 

4300  FORMAT ( 'V10X » 'CRASH, CRASH. CRASH'  ) 

WRITE (6. 4350  ) ( ROOTR ( I ) . I«1 , 4 ) 

4350  F0RnAT(/'/5X.4(E12.7.3X) ) 

GO  TO  6100 

THE  CONDENSATE  THICKNESS 

3900  DELCNI+l  )-ROOTR(l  ) 

GO  TO  4400 

4000  DELCNI+l >-R00TR(2) 

GO  TO  4400 

4100  DEL ( Nl ♦ 1 i "ROOTR ( 3 ) 

GO  TO  4400 

4200  DEL(NI+1 ) -ROOTR ( 4 ) 

4400  OEL-0  0D0 

IF(NI  NE  10)  GO  TO  4500 
URITE(6. 4600 ) 

4600  FORMAT  C ' 1 ' . ✓ /'SX . 'NP'.6X,  'X'  .12X.  ' Y' . 12X.  'T'  ) 

DO  4700  NP-l.NSNP 

4700  URITEC6.4800)  NP. X(NP  ) . V(NP ) . T(NP ) 
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ooo 


4300  FGRMAT(/2X,  13,  3X.  3CF10.6.3X  ) ) 

URIT£>6,49e0) 

4900  FORMAT < ,'/2X . ' EL ' , 3X , 'H' . 1 IX. 'EL-LENGTH' , 15X. 'Q-EL' ) 

DO  5000  KKL-l.NBOTF 
HKX-lCORCIOCL.l ) 

NICY*  ICOR  ( KKL . 2 ) 

XP*X<  MKX ) ( NKY  ) 

VP-Y<NKX)-V<NKY) 

EXY - DSQRT  < XP  **2+YPt *2 ) 

QEP- DABS  a T ( NKX  >+T ( NKY  >-2*TS  C KKL ) )*EXY*H ( KKL )/2 . OD0 ) 
QEP-QEP*DELX 

500C  WRITE (6,  5100  > KKL . Hf KKL ' . EXY . QEP 

5100  FORMAT ( '2X . 12 , 3X, E12  5 , 3X . E12  5, 10X, E12  5 ) 

WRITE' 0 5200  ) CRIT 

5200  FORMAT ( /2X, 'CONOERGENCE  CRITERIAN* ' , E15 . 8 ) 

Q FROM  THE  TOP  SIDE 

4500  DO  5300  I GEL* 1 , NEST 
KA-ICORCIGEL.l > 

KB* ICORf IGEL. 2 ) 

XQEL*X(KA)-X(KB) 

YOEL-Y(KB  )-Y(KA) 

ELM*DSGRT(XQEL#*2+YQEL**2 ) 

QEL*QEL+<  2*TS( I GEL  )-TCKA  )-T(KB  ) >*£LM*H< IGEL )/2 . 0DO 
530O  CONTINUE 

QINC(NI >*QEL*DELX 
AMTOT(NI  ) *DMTOT 
QET -QEL*DELX*NF INX2 
QT-OT+QET 
QA*QB I SOELXXNF IN*2 
QTOT-QTOT+QA 
2000  CONTINUE 

URITE  C 6 . 5400 ) HFG , NF IN . H ( NBOTF ) . TSAT . NRPM , QTOT , QT , FANGL < K I A ) 
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> 

5400  FORMAT ( ' ' .s/.SX.  'HFG-'  .E12. 5./, 5X.  'NUMBER  OF  FIM-  ' . 15. /. 5X.  'H-0* 
UT 

*-'.El2  S./.SX, 'TSAT ■ ' , E12 . 5.  ''.5X, 'RPM-' . IS./.SX, 'Q-B0T-'.E12  S./* 
.5 

*X> 'Q-TOP-'  ,E12  5./.5X. 'HALF-ANGLE- ' ,F8  3 ) 

WRITER 6 5500 ) 

5500  FORMATC '0' , 6X. 'J'.4X. 'FILM  THICKNESS' . 6X, 'Q-INCPEM' . 6X. 'MASS-TOT~ 

/ 

17X. 'TIB ' ,3X. 'TT' . 10X. 'TE' ,2X. 'TB' ) 

DO  5600  NR-l.NPIU 

5600  UPITE'6.555)  NP,DEL(NR). OB(NR ) . AMTOTfNR ) . TIB (NR ) , TTCNR  ) .TE(NR ) . T~ 
B( 

CNR ) 

555  FORMAT  ('  ' . 4X. 14, 4X.F18 . 1O.4X.F10 . 4,6X,F9  5.6X,F5 . 1 . 6X.F5 . 1 ,6X.~ 
F5 

1 . 1,6X,F5  1 ) 

WRITE f 6. 5S00  > 

5800  FORMAT ( '0' ,6X.  'J' ,6X.  'K-UALL' ,4X.  'K-FILM' ,3X,  'DENSITY' , 4X,  'UISC-~ 
FI 

CILM' -6X. 'EPSILON' ,5X, 'RADIUS', 5X. 'TBM' ,5X, 'Q-BOT' ) 

DO  5900  NG-1-HDIU.2 

5900  URITEC6 , 777  ) NG. CU(NG  ) , CF(NG  ) .RHOFCNG ) ,UF(NG ) , EPSCNG) , RING ) , TBM<~ 
NG 

* ) .QINCCNG) 

777  FORMAT  ('  ' . 4X. 14. 4X.F7 . 3. 4X,F6.4,4X,F6 . 3. 4X.F9 .7.4X.F9.7, 4X.F7  * 
5. 

15X.F5  1 . IX . 1P1D12  3 ) 

1400  CONTINUE 
6100  STOP 
END 

SUBROUTINE  COORD 
IMPLICIT  REALC8  C A-H , 0-Z ) 

COMMON/- PCRD^X ( 200  ) . Y< 200 ) , EZERO . BUIN . THICK . TALFA . APS 
COMMON/APOL/DOBF , DOTH . KFIN ( 50  > , KFF < 50  ) . IFF , JTC , JLC . JINT . KT 
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DELH-BMIM/DOBF 

xcd-o.qdq 

YU  )-THlCK:+BgiN 

H-l 

DO  321  I-l-IFF 
ZCA-KFXNCX) 
icb-kff>  i ) 

CEA-DFLOAT(ICB-ICA) 

AN-0  0D0 

DO  322  II«ICA,ICE 

X<II )«X( 1 )+N*AN*DELH*TALFA/CBA 

V(II  )-Y<l)-N*DELH 

322  AN-AN+1  ODO 
321  N-N+l 

AN-0  0D0 
ICD-ICB-ICA+1 
DO  323  J-JTC. JLC. JINT 
X(J  ) -Xf 1 > 

Y(J)-(1  0D0-AN/DOTH)*THICK 
DO  324  JJ-l.ICD 

X ( J+ J J ) »X  ( J )+ J J*EZERO.'  < 2*  ( CBA+ 1 . 0D0 ) ) 

324  Y(J+Jw)-Y(J) 

DO  32S  K-1,KT 

XC  J+JJ+K  J -X ( J+JJ  )+KXAPS/<  2 . OD0XKT ) 

325  Y(  J+-JJ+K  )"V(  J ) 

323  AN-AN+1  0D0 
RETURN 

END 

SUBROUTINE  FORHAF 
IHPLICIT  REALXS(A-H.O-Z) 

DIHENSION  B<3).C(3).EA<3.3> 

COmON/’PCRD-'X  ( 200  ) , Y ( 200  l , EZERO . BU IN . TH I CK . TALF A , APS 
connoN/riAFO/A  ( 200 . so ) . f ( 200 . 1 ) . h c 200 ) . ts  < 200  > . tsat  . ck  . nel  . nsnp  . n* 

BA 
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*N, I COR <200. 3) 

DO  101  N-l.NSNP 
F(N. 1 i-0  ODO 
DO  102  NA-l.NBAN 

102  A<N.  1*1(4  )-0 . 0D0 

101  CONTINUE 

DO  103  IEL-1 . NEL 
IA-ICOR' IEL, 1 ) 

IB-ICOR*.  IEL- 2) 

IC-ICOR( IEL. 3 ) 

B(  1 IB  >-V(  10  ) 

E<2  )■*.’(  IC  )-Y(  IA  ) 

Bt3 )-Vl Irt )-Y< IB  ) 

C( 1 )-X( IC  >-X( IB  ) 

C<2>-XlIA)-X(IC » 

C(3  )■><(  IB  )-X(  IA  ) 

EL-DSGPT(C<3  )X*2+B<  3)**2> 

AS-DttBS< ( B ( 1 )*C(2)-B(2)*C(1 ) )/2.0D0) 
HC-H(  IED/CK 
DO  104  J-1.3 
JJ-ICORC IEL. J ) 

DO  105  K-1.3 
KK-ICOR(IEL.K) 

EA<J,K)-(B( J)XB(K>+C<J)*C<K))/(4XAS> 

IFCHC  EQ  0 0D0 ) GO  TO  106 

HEL-HCXEL/6  0D0 

IF( J EQ  3)  GO  TO  106 

IF(< . EQ  3)  GO  TO  106 

IFCJ.EQ  K>  GO  TO  1C7 

EA( J,K  )*EA( J.K  )+HEL 

GO  TO  106 

107  EA<J.K)«EA<J.K)+2*HEL 

106  IFCKK  LT.JJ)  GO  TO  105 
NU-KK-JJ+1 
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A < J J , NU  ' •*  A < J J . Nl J ) + E A < J . K. > 

105  CONTINUE 
104  CONTINUE 

FE-HC*TS(IEL>*EL'2  ODO 
F<  IA.  1 *F<!  IA • 1 >+FE 
F ( ID/ 1 '*F ( IE . 1 ’+FE 
103  CONTINUE 
RETUPN 
END 

SUBROUTINE  EANPEC ( NEQ • NAXB . NUtC ) 

IP1PLICIT  REAL*3(A-H.0-Z' 

COMMON  'PCRD'  rtf  200  > , V(  20o  ) . E2ER0 . BUIN . THICK . TALFA . APS 
COnnON/flrtFO/ A ( 200 . 50  ) . F < 200  - 1 > . H ( 200  ) , TS  < 200 ) . TSAT . CK . NEL . NSNP . N~ 
BA 

*N, IC0Rf2«0.3> 

LOOP-NEQ-1 
DO  201  1*1. LOOP 
riB-i+i 

NB-NINOf I+NAXB-1 , NEQ  ) 

DO  201  J-MB.NB 

L-J+2-MB 

D-A(I.L)/A(I.l) 

• DO  202  nri-l.NVEC 

202  F(j.pin)-F(j.m)-D*F(i.nn) 
m-niNOcriAXB-Lt-i  ,neq-j+i  ) 
do  201  K-i.nn 

NN-L+K-1 

201  A(J.K)-A(J.K)-D*A(I.NN) 

DO  203  I-l.NUEC 

203  F(NEQ. I )-F(NEQ. I )/A(NEQ.l ) 

DO  204  I -2. NEQ 
J-NEQ-I+1 

K-NINO  C NEQ- J+l . P1AXB  ) 

DO  204  m-l/NUEC 
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205  Fi'J  NM  >«F  ( J . MM  '-A*.  J • L >tF ‘.MB , MM ) 


20-4  F<  J,NN  '-F(  J.MM  J^'ACJ.  1 ) 
RETURN 
EMD 
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